Publisher Correction: Mapping propagation of collective modes in Bi2Se3 and Bi2Te2.2Se0.8 topological insulators by nearfield terahertz nanoscopy

Mapping propagation of collective modes in Bi2Se3 and Bi2Te2.2Se0.8 topological insulators by near-field terahertz nanoscopy Eva Arianna Aurelia Pogna, Leonardo Viti, Antonio Politano, Massimo Brambilla, Gaetano Scamarcio, Miriam Serena Vitiello NEST, CNR-Istituto Nanoscienze and Scuola Normale Superiore, Piazza San Silvestro 12, 56127 Pisa, Italy Department of Physical and Chemical Sciences, University of L’Aquila, via Vetoio 1, 67100 L’Aquila, Italy Dipartimento Interateneo di Fisica, Università degli Studi e Politecnico di Bari, via Amendola 173, 70126 Bari, Italy


Supplementary Note 1: Micro-Raman spectroscopy
In order to evaluate the crystalline quality of the exfoliated topological insulator (TI) Bi 2 Se 3 and

Spatial resolution and far-field background
In our nanoscopy experiments, the tip is operated in tapping mode with oscillation frequency Ω.
The near-field signal depends non-linearly on the tip-sample distance as it is shown in Supplementary Fig. 3.
The near-field interaction decreases over distances comparable to the tip radius (tens of nanometers), whereas the background varies on the radiation wavelength scale (tens of microns for THz frequency light). Accordingly, for tapping amplitudes < 500 nm, the background slowly changes and it contributes mostly to the dc term and to the lower harmonics.

Supplementary Figure 3: Approach curves for hyperspectral nanoimaging measurements. a-b
Approach curves at the harmonics of the tapping frequency of order n= 2, 3, 4 recorded with TDS system on gold for fixed time delay corresponding to an amplitude peak of the THz pulse, reported in semilogarithmic scale (a); in linear scale (b) to calculate the distance h 1/2 , at which the normalized signal drops to half of its maximum. Inset: h 1/2 plotted as a function of the order n of the harmonic.
To find a good compromise between signal intensity and suppression of the background due to non-local scatterers, we select the harmonic n= 2 to analyze the near-field maps, being the lowest order Fourier component at which the signal reaches the noise level once retracting the tip from the sample in Supplementary Fig. 3. At n= 2 and higher order harmonics, the signal reaches the same value for large tip-sample distance, corresponding to the noise level of our detection system.
The decay distance h 1/2 , at which the signal amplitude decays by a factor of 2, is equal to 63 nm at n= 2 and decreases with the demodulation order n, due to the virtual tip sharpening effect 3 . We then characterize the signal amplitude s n and the lateral resolution Δx by analyzing the line profiles across the interface between an Au marker and the Si substrate at a fixed time delay (see Supplementary Fig. 4).

Supplementary Figure 4: Spatial resolution of hyperspectral nanoimaging measurements. a
Topography scan of Au marker along a line orthogonal to its edge, taking the substrate as z = 0 reference. b THz near-field self-mixing signal at the second harmonic s 2 (black dots) from the hyperspectral nanoimaging experiment evaluated at a fixed time delay on the same line as panel (a) together with the best fit using sigmoidal function (solid blue line). The first derivative of the signal ds 2 /dx (black dots) is reported in the inset, together with Lorentzian fit function (blue line) used to estimate the lateral spatial resolution Δx.
We define the spatial resolution, in analogy to other microscopy techniques, as the full width at

Multilayer scattering model
In order to disentangle the dielectric response from the bulk and the surface, we combine  Fig. 6c, 6e) Detectorless self-mixing interferometric nanoscopy is performed coupling our THz-QCL sources to a commercial s-SNOM from Neaspec/attocube. The spatial resolution of the self-mixing interferometry experiment is evaluated through the analysis of the signal variation at the edge of a 50 nm high Au marker (Supplementary Fig. 8).

Self-mixing interferometry near-field maps
The first derivative of the step-like line profile in Fig. S8b is then fitted with a Lorentzian function with FWHM of 30 ± 5 nm.

Self-mixing fringes analysis
Supplementary  Fig. 9d-f). The position of the edge corresponds to x= 0, for x< 0 the tip is on the substrate and for x> 0 on the flake.
The phase jump is almost thickness independent in Bi 2 Se 3 , while it shows a visible dependence from the flake thickness in Bi 2 Te 2.2 Se 0.8 only at the highest pumping frequency (89.7 cm -1 ). In all cases, although the retrieved trends do not perfectly match those expected (Fig. 1d, main text) For both materials, we use these trends to evaluate the phase variation Δφ 3 between substrate and flake reported in Fig.3 of the main text. Error bars are the 95% confidence bands from the sinusoidal fit of the self-mixing fringes.

Signal at the flake edges
Supplementary Figures 10 and 11 show the n= 3 near-field signals collected close to the edge of flakes of Bi 2 Se 3 ( Supplementary Fig. 10) and Bi 2 Te 2.2 Se 0.8 ( Supplementary Fig. 11) of different thickness, with QCL sources emitting at 76.7 cm -1 and 89.7 cm -1 . The near-field self-mixing maps are analyzed to identify fingerprints of interference patterns due to the superposition of propagating polariton modes launched by the tip and reflected at the flake edges. Line profiles are extracted from the maps along cuts with direction orthogonal to the edges. The signal oscillations are analyzed with the function s 3 described in the main text to extract the periodicity λ p shown in Figure 5, as a function of the flake thickness and of the laser frequency. To fit the data, we exclude the first highest peak at the flake edge, which can include contributions arising from the edge modes. The oscillations appear damped (only less than 3 periods can be observed) and they show higher amplitude for a signal demodulation order n= 3 than for higher orders.

Bulk and Dirac Plasmon dispersion in Bi 2 Te 2.2 Se 0.8
Supplementary Fig. 12 shows the simulated dispersion for bulk and Dirac plasmons in Bi 2 Te 2.2 Se 0.8, calculated following the same procedure described in the main text for Bi 2 Se 3 , and employing, as physical parameters, the ones commonly adopted for Bi 2 Te 3 . Specifically, we consider an effective mass 11 m eff = 0.044 m 0 , where m 0 is the electron mass in vacuum, and a Fermi velocity v=0.5×10 8 cm/s, which are different from those adopted in Bi 2 Se 3 (m eff = 0.15 m 0 , v=0.623 ×10 8 cm/s 14 ). The simulated q(λ p ) reported in Supplementary Fig. 12 are at least one order of magnitude smaller (larger) than those observed in experiments.
For calculating the dispersion of bulk plasmons, we consider the bulk carrier density previously evaluated for similar stoichiometry 11 n bv = 3×10 18 cm -3 . The total carrier density n M is then evaluated by considering a square-root dependence on the flake thickness d as n M = n s + C⋅n bv ⋅ d 1/2 (Ref. 15). Sheet carrier density is taken as n s = 2×10 12 cm -2 , estimated by multiplying n bv from Ref.
11 for the thickness of the flake used for the n bv determination (6 QL~8nm). Figure 12: Plasmon dispersion. a-b Predicted plasmon energy dispersions in Bi 2 Te 2.2 Se 0.8 , calculated using the physical parameters measured experimentally 9 in Bi 2 Te 2.4 Se 0.6 , for flakes of various thickness ranging from 118 nm to 19 nm for massive bulk plasmons (a) and Dirac plasmons (b); c-d predicted plasmon wavelength λ p = 2π/q calculated using the q values of panels (a-b), for massive bulk plasmons (c) and Dirac plasmons (d).

Bulk dielectric permittivities of Bi 2 Se 3 and Bi 2 Te 3
The two material systems investigated in the present work are modeled as anisotropic Drude-Lorentz materials with two different bulk permittivity functions for the in-plane directions ε ! ω , parallel to the basal plane, and the out-of-plane direction ε // ω , along the c-axis. The giant anisotropy of the crystal structure and of the vibrational properties gives rise to a strong anisotropy in their dielectric response of the two investigated TIs. Each of the two permittivity components is here described by two oscillators corresponding to the α-and β-phonons, while the electronic contribution in the frequency range of interest, which is well-below the electronic gap (ω < E g /ħ ~2400 cm -1 in Bi 2 S 3 ), is described by ε ! .
In Supplementary Table 1 and Table 2 we report the parameters used to model Bi 2 Se 3 and Bi 2 Te 3 respectively, chosen from available experimental and theoretical literature [16][17][18] .  The near-field scattering signal expected for bulk Bi 2 Se 3 and Bi 2 Te 3 reported in Supplementary   Fig.5a has been obtained considering the bulk permittivity in Supplementary Fig. 13 and applying the model for the near-field interaction 17 .

Simulation of the near-field signal s 3
We model the tip as a metallic spheroid with curvature radius r= 20 nm and length l= 80 µm considering the scattering problem in the quasi-static approximation 16,18 . We compute the full scattering signal s for a set of different tip heights z tip , defined with respect to the sample surface.
The tapping is modelled as an oscillation of z tip with amplitude Δz= 160 nm with minimum tipsample distance z min = 10 nm. We calculate s as the product of the tip-sample polarizability χ ⊥ (ω, z tip ) and the far-field factor F= |1+r p | 2 . The far-field factor accounts for the fact that the source and the detector are in the far-field [18][19][20] and includes the reflectivity r p evaluated at q= c/ω. The thirdorder near-field signal s 3 corresponds to the third Fourier harmonic of s(z tip ). The third order polarizability χ 3 ⊥ (ω, z tip ) depends both on the geometrical properties of the tip and on the sample reflectivity r p , and, as a consequence of the finite q-dependence of r p of the TI, it is momentum dependent. . Specifically, we compute the k residues R k ⊥ and poles β k ⊥ of the polarizability χ 3 as a function of z tip , following the generalized spectral method for spheroidal tips 18 including k= 200 terms, and using an effective reflectivity averaged within the momentum range from 0 up to 3 nm -1 , weighting the different q with normalized Gaussian distribution functions centered at 1/z tip. The retrieved signal for a slab of Bi 2 Te 3 with chemical potential µ= 200 meV (i.e. below the bandgap of Bi 2 Te 2.2 Se 0.8 ) is reported in Supplementary Fig.   15 after dividing for the signal expected from the SiO 2 /Si substrate modelled by constant permittivity ε sub = 6.3. We observe a peak in s 3 at the frequency of the phonon mode ω ! ! = 50 cm -1 (see Supplementary Note 4.5 for bulk phonon modelling) which is the analogue of the s 3 peak at 64 cm -1 previously predicted for Bi 2 Se 3 16 . This peak originates from the far-field factor that has a maximum at ω ! ! . Conversely, the peak centered at 100 cm -1 can be attributed to surface modes: hyperbolic phonon-polaritons that hybridize with Dirac plasmons. The coupling with the tip, included in the tip-sample polarizability, results into a small (few cm -1 ) difference between the frequency position of this peak in s 3 and the corresponding peak in the sample reflectivity Im (r p ) at ω α // = 95 cm -1 . The simulated trend accounts for the contrast change in Fig. 3d when passing from ω= 66.7 cm -1 to 76.7 cm -1 , and predicts the signal increase observed at 89.7 cm -1 as approaching the expected peak centered at 100 cm -1 and the stronger thickness dependence at 66.7 cm -1 compared to 89.7 cm -1 that results from the decrease of r p with thickness. Finally, the signal phase is expected to increase approaching the peak in the amplitude at 100 cm -1 and could be compatible with the phase variation observed at 90 cm -1 in the experiments, see Supplementary   Fig. 9, after considering a rigid shift to lower frequency of the predicted trend, shown in Fig.   Supplementary Fig. 15b.